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I. INTRODUCTION 

Solid state spin nanodevices are known as very promis- 
ing candidates for quantum computatio n 1 ' 2 "^ 4 and also 
for quantum communication 5 '. They provide a scalable 
system that can easily be integrated into standard sili- 
con technology. A drawback of such systems compared 
to other proposals for quantum computing, such as ion 
traps& and cavity QEDi, are the many degrees of free- 
dom of the surrounding material causing dissipation and 
decoherence^. The first step in overcoming these disad- 
vantages is to be able to model them. 

An important contribution to quantum noise in solid 
state systems arises from the nuclear spins, and recently 
much work has been devoted to the modeling of spin 
bath system ofti 1 fti 11 i 12 i 1 i?i 14 i 1 ''?i 1 fii 1 T (for a review, see Refs. 

The interaction of a central spin with a bath of 
environmental spins often leads to strong non-Markovian 
behavior. The usual derivations of Markovian quantum 
master equations known, e. g., from atomic physics*^ and 
quantum optics^ therefore fail for many spin bath mod- 
els and a detailed investigation of methods is required 
which are capable of going beyond the Markovian ap- 
proximation. 

In this paper, we examine the reduced dynamics of a 
simple spin star system. The advantage of this model 
is that, while showing several interesting features such 
as partial decoherence and strong non-Markovian be- 
havior, it is exactly solvable due to its high symme- 
try. The model therefore represents an appropriate ex- 
ample for a general discussion of the performance of 
various non-Markovian methods. We study here the 
Nakajima-Zwanzig^ 2 * 2 ^* 2 ^ and the time-convolutionless 
projection operator technique 2 * 5 *2£* 2 i 2 £ and derive and 
analyze the perturbation expansions of the corresponding 
non-Markovian master equations. 

The paper is organized as follows. In Sec.|Ul we intro- 
duce the model investigated, a spin star model involving 
a Heisenberg XX coupling (Sec. Ill All , and determine the 
exact time evolution of the central spin fSec. Ill 5)) , In 



Sec. Ill CI we analyze the limit of an infinite number of 
bath spins, discuss the behavior of the von Neumann 
entropy of the central spin, and demonstrate that the 
model exhibits complete relaxation and partial decoher- 
ence. Several non-Markovian approximation techniques 
are discussed in Sec. lIIII The dynamic equations found in 
second order in the coupling are introduced in Sec. lIII Al 
where it is also demonstrated that the prominent Born- 
Markov approximation is not applicable to the spin star 
model. Employing a technique which enables the cal- 
culation of the correlation functions of the spin bath, 
we derive in Sec. IIIIBl the perturbation expansions cor- 
responding to the Nakajima-Zwanzig and to the time- 
convolutionless projection operator technique and com- 
pare these approximations with the analytical solution 
for the dynamics of the central spin. Finally, the conclu- 
sions are drawn in Sec. lIVI 



II. EXACT DYNAMICS 
A. The model 

We consider a "spin star" configuration^ 9 , which con- 
sists of N + 1 localized spin-i particles. One of the spins 
is located at the center of the star, while the other spins, 
labeled by an index i = 1, 2, . . . , TV, surround the central 
spin at equal distances on a sphere. In the language of 
open quantum systems*^ we regard the central spin with 
Pauli spin operator <r as an open system living in a two- 
dimensional Hilbert space Hs and the surrounding spins 
described by the spin operators cr^ as a spin bath with 
Hilbert space Hb which is given by an iV-fold tensor 
product of two dimensional spaces. 

The central spin er interacts with the bath spins erM 
via a Heisenberg XX interaction^ represented through 
the Hamiltonian 

aH = 2a(er+J_ +cr_J+), (1) 
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and 
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represents the raising and lowering operators of the zth 
bath spin. The Heisenberg XX coupling has been found 
to be an effective Hamiltonian for the interaction of some 
quantum dot systems^!. Equation describes a very 
simple time independent interaction with equal coupling 
strength a for all bath spins. It is invariant under rota- 
tions around the z-axis. The operator J = \ J2iLi a ^ 
represents the total spin angular momentum of the bath 
(units are chosen such that h = 1). The central spin thus 
couples to the collective bath angular momentum. 

We introduce an ON basis in the bath Hilbert space 
Hb consisting of states \j, to, v) . These states are defined 
as eigenstates of J3 (eigenvalue to) and of J 2 [eigenvalue 
j{j + 1) ]• The index v labels the different eigenstates in 
the eigenspace M.j, m belonging to a given pair (j, to) of 
quantum numbers. As usual, j < y and —j<m<j. 
The dimension of M.j. m is given by the expression 29,32 



n{j,N) 



N 
N/2 - j 



N 

N/2-3-I 



(4) 



We further introduce the usual superoperator notation 
for the Liouville operator 



Cp(t) = -i[H,p(t)] 



(5) 



where pit) denotes the density matrix of the total system 
in the Hilbert space Tis <S> Tis- The formal solution of 
the von Neumann equation 



— p(t) = aCp(t) 



can then be written as 



p(t) = exp(aCt) p(0). 



(6) 



(7) 



Our main goal is to derive the dynamics of the reduced 
density matrix 



p s (t) = tT B {p(t)} , 



(8) 



where trs denotes the partial trace taken over the Hilbert 
space Kb of the spin bath. The reduced density matrix 
is completely determined in terms of the Bloch vector 

«i(f) 

v(t) = I v 2 (t) I =tTs{<rpsffl (9) 

«s(*) 



through the relationship 

ps (t)- 1 ( 1 + Vs W ~ iv2 ^ 



2 \ vi(t) + iv 2 (t) l-«a(t) 



, (io) 



where trg denotes the partial trace over the Hilbert 
space Tis of the central spin. We note that the length 
r(t) = \v(t)\ of the Bloch vector is equal to 1 if and only 
if ps(t) describes a pure state, and that the von Neu- 
mann entropy S of the central spin can be expressed as 
a function of the length r(t) of the Bloch vector: 



S = trs {-ps^nps} 



(11) 



= In 2 - -(1 - r)ln(l - r) + -(1 + r) ln(l + r). 

The initial state of the reduced system at t = is taken 
to be an arbitrary (possibly mixed) state 



Ps(0) 



i±fM v _(o) 
v + (0) ^ 



(12) 



while the spin bath is assumed to be in an unpolarized 
infinite temperature state: 



Pfl(0) 



2- N I F 



(13) 



Here, Is denotes the unit matrix in Hb, and we have 
defined the v± as linear combinations of the components 
vi 2 of the Bloch vector: 



v± 



(14) 



The initial state of the total system is given by an uncor- 
rected product state /?s(0) <g> /5_b(0). 



B. Reduced System Dynamics 

In this section, we will derive the exact dynamics of the 
reduced density matrix ps(t) for the model given above. 
One possibility of obtaining the evolution of the central 
spin is to substitute Eq. @ into Eq. JHl and to expand 
the exponential with respect to the coupling. This yields 

Ps it) = tr B {exp(a£t)p s (0)®PB(0)} 

= £^-^{£^(0)8 2-^}. (15) 



It is easy to verify that we have 

H 2n = 4" [<j + a-(J-J + ) n + a_a + (J+J_) n ] (16) 

and 

H 2n+i = 2 . 4 « [ a+( T_(J + J_) n + a-a + {J-J + ) n ] . (17) 

We note that such simple expressions are obtained since 
a term 03 J3 is missing in the interaction Hamiltonian. 
We substitute the last two equations into 

£ n p = i n jr(-l) k ( n \H k P H n - k (18) 

k=Q ^ ' 
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(19) 



to get the formulas 

ti B {C 2k - 1 ps(0)(g>2- N I B } =0 
and 

tru {£ 2k Ps (0) ® 2- N I B ) = (-16)" 0^(73 ( 20) 



si=0 



which hold for all k = 1, 2, — Here, we have introduced 
the bath correlation functions 



1? 



fe 
fc-i 



^tr s {(J+J_) fc }, (21) 
^tr B {(J + J_) fe -'(J_J+)'}. (22) 



We will come back to these correlation functions when 
we discuss approximation techniques in Sec. IIIII 

Using the formulas ltl9)l and l(20l) in Eq. ltl5)l we can 
express the components of the Bloch vector as follows, 



V±(t) = /l2(t)«±(0), 

v 3 (t) = f 3 (t)v 3 (0), 
where we have introduced the functions 

/l2(t) = 

tr B | cos 2 ^/ J + J_ at cos 2\f~J~J+cd 

and 

/3(f) = trs I cos 4-\/ J+J-at 



(23) 
(24) 



(25) 



'4 



(26) 



Calculating the traces over the spin bath in the eigenbasis 
of J3 and J 2 using 



Jt J ± \h m , v ) = U Tm)(j±m + 1) \j, m, v) 
we find 

Mt) = 



(27) 



cos [2/i (j, m)at] cos [2/i(j, — m)crf] 
2^ ' 



3,m 



and 



(28) 



(29) 



3,m 



where we have introduced the quantity h(j, m) = 

Vti + m )U -m+1). 

Thus we have determined the exact dynamics of the 
reduced system: The density matrix ps(t) of the central 
spin is given through the components of the Bloch vector 
which are provided by the relations l12.'ill . H24I1 and l|28[). 



II29H . We note that the dynamics can be expressed com- 
pletely through only two real- valued functions /12 (t) and 
fz{t). This fact is connected to the rotational symmetry 
of the system. 

The reduced system dynamics has been obtained above 
with the help of an expansion of exp(a£t) with respect 
to the coupling constant a. It should be clear that, alter- 
natively, the behavior of the central spin may be found 
directly from the solution of the Schrodinger equation 
for the composite system. This solution can easily be 
constructed by making use of the fact that the subspaces 
spanned by the states |+)<8>|j, m, v) and |— m + 1, v) 
are invariant under the time evolution, where |±) denotes 
the eigenstate of (T3 belonging to the eigenvalue ±1. 

Sometimes it is useful to express the reduced dynam- 
ics in terms of superoperators instead of the Bloch vec- 
tor. To this end, we introduce superoperators S± and S3, 
which are defined by their action on an arbitrary operator 
A: 

S ± A = a±Aa T - - {a T <i ± ,A} , (30) 



S 3 A ee a 3 Aa 3 - A. 



(31) 



With these definitions we may write the reduced density 
matrix as follows, 

Is 
2 



Ps{t) 



1 



/3W-/12W S3-/3W (S++S-) 



(32) 
Ps(0), 



where Is denotes the unit matrix in Tls- Due to the 
non-unitary behavior of the reduced system, the super- 
operator on the right-hand side is not invertible for all 
times. This point will become important later on when 
we discuss approximation strategies. 



C. The Limit of an infinite number of bath spins 

The explicit solution constructed in the previous sec- 
tion takes on a relatively simple form in the limit N — > 00 
of an infinite number of bath spins. It is demonstrated 
in the appendix that for large TV the bath correlation 
functions approach the asymptotic expression 



R 



k-l 



2 k 



(33) 



Consequently, in Eq. 11511 a term of the order N k always 
occurs together with a factor of a 2k . A non-trivial finite 
limit N — > 00 therefore exists if we rescale the coupling 
constant as follows, 



a 



N 



(34) 



Using this approximation in Eq. 112011 one obtains 
ti- B {C 2k p s (0)® 2- N I B ) 

« ( ~ 8j p fc! M0)<7 3 + MOW + «+(0)tr_) . (35) 



4 



If we insert this into Eq. 11 Till we get an infinite series 
which yields the following expressions for the functions 

/ia(t) and f 3 (t), 



f 12 (t) = l+g(t), f 3 (t) = l + 2g(t), (36) 



where 



g(t) = -aiexp(-2a 2 i 2 )y^erfi(\/2crf). (37) 

Note that erfi(cc) is the imaginary error function. It is a 
real-valued function defined by 



erfi(x) 



ert(ix) 



which leads to the Taylor series 



„2fc+l 



V^fc!(2fc + 1) 



(38) 



(39) 



Figures and El show that this approximation obtained 
in the limit of an infinite number of bath spins is already 
reasonable for TV w 200. 



v 3 (t) 



v 3 (t) 



20 bath spins N -> oo 



exact 




5 ,--2 "23 "3 
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v ± (t) 
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Figure 1: Comparison of the limit TV — > oo [see Eqs. i'M)l and 
G3J] with the exact functions for TV = 20 and TV = 200. The 
plot shows the f±-component of the Bloch vector. 



By contrast to erf (a;), the imaginary error function is 
not bounded. However, g(t) is bounded, and in the limit 
t — > oo we have g(t) — > — 4. Thus, in this limit the system 



Figure 2: Comparison of the limit TV — > oo [see Eqs. ft-itijl and 
JS}] with the exact functions for TV = 20 and TV = 200. The 
figure shows the «3-component of the Bloch vector. 



is described by the stationary density matrix 



lim lim ps(t) 

t^oo N— too 




(40) 



The 3-component of the Bloch vector relaxes to zero, 
while the off-diagonal elements of the density matrix 
show partial decoherence, i.e. they assume half of their 
original values. This behavior is also reflected in the von 
Neumann entropy of the reduced system. Its dynam- 
ics depends on the initial entropy parametrized by r(0) 
and on t>3(0), which is obvious because the entropy is a 
scalar quantity and the system is invariant under rota- 
tions around the z-axis. Figure shows the entropy as a 
function of time for different initial conditions. 



III. APPROXIMATION TECHNIQUES 

In this section we will apply different approximation 
techniques to the spin star model introduced and dis- 
cussed in the previous section. Due to the simplicity of 
this model we can not only integrate exactly the reduced 
system dynamics, but also construct explicitly the vari- 
ous master equations for the density matrix of the central 
spin and analyze and compare their perturbation expan- 
sions. In the following discussion we will stick to the 
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S(t) 




r(0) = 0.7 
v 3 (0) = 0.3 



Figure 3: Von Neumann entropy S(t) of the reduced sys- 
tem for different initial conditions in the limit N —* go. 
Smax = In 2 is the maximal entropy for a qubit, represent- 
ing a completely mixed state. 



Bloch vector notation. Each of the master equations ob- 
tained can easily be transformed into an equation involv- 
ing Lindblad superoperators [see Eqs. l(30)) . (1311) ] using 
the translation rules 

V303 = I -5 3 - 5+ - 5_ Us, (41) 



V+<7_ + V-(J-\ 



-SzPs- 



(42) 



A. Second order approximations 



The second order approximation of the master equa- 
tion for the reduced system is usually obtained within the 
Born approximations,. It is equivalent to the second order 
of the Nakajima-Zwanzig projection operator technique, 
which will be discussed systematically in Sec. ITTTBl In 
our model the Born approximation leads to the master 
equation 



hit) 



I dstT B {[H,[H,p s (s)® p B 
Jo 



(0)]]} 



-8c^Qi / ds (v+(s)a- + v-(s)cr+ + v 3 (s)a 3 ) , (43) 
'0 



where the bath correlation function is found to be 
1 

1 



Qi 



-tT B {J+J-} 



2 N 



tr B 



It is important to notice that Qi, as well as all other 
bath correlation functions are independent of time. This 
is to be contrasted to those situation in which the bath 
correlation functions decay rapidly and which therefore 
allow the derivation of a Markovian master equation. 
The time-independence of the correlation functions is the 
main reason for the non-Markovian behavior of the spin 
bath model. 

The integro-differential equation 114311 can easily be 
solved by a Laplace transformation with the solution 

v±(t) 



«±(0) 
v 3 (0) 



cos 



foVNat) 



(45) 



(46) 



In many physical applications the integration of the 
integro-differential equation is much more complicated 
and one tries to approximate the dynamics through a 
master equation which is local in time. To this end, the 
terms v±(s) and v 3 (s) under the integral in Eq. 14311 are 
replaced by v±(t) and v 3 (t), respectively. We thus arrive 
at the time-local master equation 



dt 



Ps(t) 

= -4Na 2 I ds (v + (t)a- + V-(t)a + + v 3 (t)a 3 ) 
Jo 

= -4Nta 2 (v+(t)a- + V-(t)a+ + v 3 {t)<r 3 ) , (47) 

which is sometimes referred to as Redfield equation. Also 
this master equation is easily solved to give the expres- 
sions 



v±(t) 
«±(0) 

«3(0) 



exp(-2Na 2 t 2 ), 
exp(-4iVa 2 i 2 ). 



(48) 
(49) 



The Redfield equation is equivalent to the second order of 
the time-convolutionless projection operator technique, 
which will also be discussed in detail in Sec. lIIIBl 

In order to obtain, finally, a Markovian master equa- 
tion, i.e. a time-local equation involving a time indepen- 
dent generator, one pushes the upper limit of the integral 
in Eq. 114 7J1 to infinity, that is one studies the limit t — > 00 
of the master equation. This limit leads to the Born- 
Markov approximation of the reduced dynamics. In the 
present model, however, it is not possible to perform this 
approximation because the integrand does not vanish for 
large t. Thus, the Born-Markov limit does not exist for 
the spin bath model investigated here and the descrip- 
tion of relaxation and decoherence processes requires the 
usage of non-Markovian methods. 



1 

2N 
N 



tr B 



(44) 



B. Higher order approximations 

A systematic approach to obtain approximate non- 
Markovian master equation in any desired order is pro- 
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vided by the projection operator techniques. We define 
a projection superoperator V through the relation 



and 



Vp = tr B {p} <8> Pb 



(50) 



with the reference state ps — Pb{0) and introduce the 
notation 



(X) = vxv 



(51) 



for any superoperator X . Note that the "moments" (X n ) 
are operators in the total Hilbert space lis ® Ti-B of the 
combined system. 

There are two main projection operator methods: 
The Nakajima-Zwanzig (NZ) technique and the time- 
convolutionless (TCL) technique. In our model, the ini- 
tial conditions factorize. The NZ and the TCL method 
therefore lead to relatively simple, homogeneous equa- 
tions of motion. The NZ master equation is an integro- 
differential equation for the reduced density matrix with 
a memory Af(t, r), which takes the form 



Ps{t)®p B = drM{t,T)p s {T)® p B , (52) 
Jo 

while the TCL master equation is a time-local equation 
of motion with a time-dependent generator K(t), which 
reads 



ps(t) <8> PB= rC(t)ps(t) <g> PB- 



(53) 



Both the NZ and the TCL master equation can of course 
be expanded with respect to the coupling strength a. 
Since the interaction Hamiltonian is time independent, 
this expansion yields 

/ drAf(t,r)p s (r) = J2^ n Mt,r){C n ) pc ps(r) (54) 

J ° n=l 



for the NZ master equation, and 



\ (n-i): 



(£«) 



(55) 



for the TCL master equation, where we have introduced 
the integral operator 

T„(t, r) = dh r 41 dt 2 ■ ■ ■ J*- 3 dtn-2 ft' 2 dT ( 56 ) 

for the ease of notion. The symbol {C n ) pc denotes the 
partial cumulants and (£ n ) oc the ordered cumulants of 
order n. Their definitions can be found in Refs. l25Tl2?l 
In our model we have 



(C 2n+1 ) =(£ 2n+1 ) 

I pc \ / c 







<£ 2 > 

\ 1 pc 


= <£ 2 >, 






(c 2 ) 

\ 1 oc 


= <£ 2 >- 






(C A ) 

\ 1 pc 


= (O- 






\ 1 OC 


= <£ 4 >- 


-3<£ 2 ) 2 , 




\ 1 pc 


= 


-2<£ 2 )(£ 4 )4 


-3<£ 2 ) 3 , 


(c 6 ) 

\ 1 oc 


= 


-15(£ 2 )<£ 4 ) 


+ 30(£ 2 ) 3 



In the time independent case the ordered cumulants are 
just the ordinary cumulants know from classical statis- 
tics. To calculate theses functions one can again use Eq. 
l!2(H . The functions Qfc and Rf are real polynomials in 
TV of order k. A method of determining these polynomi- 
als is sketched in the appendix. 

If we express the resulting master equations in terms 
of v±{t) and v^{t), we get for the TCL technique 



TCL: v±(t) 



2n s ^nt 



2n-l 



\n— 1 v ' / 



2n ^Q2nt 



2n-l 



J2* 2n =Z^r)^(t), (58) 



(2n-l)! 
and for the NZ method 

NZ: v±(t) = ^a 2n S 2 n2"n(i,T)^ ± (r), (59) 

«s(t) = ^a 2n 2g 2n 2- n (t,r)^ v 3 (r). (60) 

The quantities S2n, «2n, ?2n and qi n represent real poly- 
nomials in TV of the order n. For example, we have 

q 2 = -47V, 
q 4 = -32TV 2 , 

q 6 = -1024TV+ 1536TV 2 - 1536TV 3 , 

q 2 = -47V, 
q\ = 32TV 2 , 

q 6 = -10247V + 1536TV 2 - 12807V 3 , 

s 2 = -47V, 

s 4 = -487V + 167V 2 , 

s 6 = -10247V- 3847V 2 + 3847V 3 , 



s 2 = -47V, 

s 4 = -487V + 487V 2 , 

s 6 = -10247V + 21127V 2 -12167V 3 , 
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The (2n)th-order approximation of the master equa- 
tions (denoted by TCL2n and NZ2n, respectively) is ob- 
tained by truncating the sums in Eqs. llfiTjl . ll58|) and in 
Eqs. Ipl after the nth term. In the TCL case, 

the resulting ordinary differential equations can be inte- 
grated very easily. The equation of motion of the NZ 
method can be solved with the help of a Laplace trans- 
formation. However, it may be very involved to carry out 
the inverse transformation for higher orders. For exam- 
ple, the solution of the twelfth order of the NZ equation 
as obtained by standard computer algebra tools is filling 
some hundred pages, whereas the solution of the TCL 
equation can be written in a single line. 

The solutions of the master equations in second and 
fourth order are plotted in Figs. and El together with 
the exact solutions. We observe that both methods lead 
to a good approximation of the short-time behavior of 
the components of the Bloch vector. We further see that 
the TCL technique is not only easier to solve, but also 
provides a better approximation of the dynamics within 
a given order. 




0.05 \ 0.1 0.15 0.2 



at 




0.05 0.1 0.15 0.2 



at 



-0.5 
-1 



Figure 4: Comparison of the TCL and the NZ technique with 
the exact solution. The plot shows the approximations to 
second and fourth order in a and the exact solution of v± (t) 
[see Eqs. i'2'31 and 1281 1 for a bath of 100 spins. 



Since the TCL and the NZ method lead to expansions 
of the equations of motion and not of their solutions, the 
solutions of the truncated equations may contain terms 
of arbitrary order in the coupling strength. For example, 
even though TCL2 is a second order approximation, the 




-0.5 
-1 



Figure 5: Comparison of the TCL and the NZ technique with 
the exact solution. The plot shows the second and the fourth 
order approximations as well as the exact solution of V3(t) 
[see Eqs. 124II and l|29ll ] for a bath of 100 spins. 



corresponding solution given by Eqs. 148|l and jHH) con- 
tains infinitely many orders. Of course, the expansion of 
the exact solution coincides with the expansion of the ap- 
proximations obtained with TCL2n or NZ2n within the 
(2n)th order. The error of TCL2 or NZ2, for example, is 
therefore a term of order a 4 , as is illustrated in Fig. El 



0.01 
0.0001 
1.x10 -6 
1.x10" 8 



E(t) 



TCL 2 



NZ 2 /,■ 



0.001 



0.01 



ot 



Figure 6: Error of TCL2 and NZ2: The plot shows the de- 



viation e(£) = «±(£) — v\ 



*(t)| of the exact solution v±(t) 



from the approximate solution w^ pprox (i) for small at. 



Concerning the long-time behavior, both the TCL and 
the NZ method may lead to very bad approximations. 
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For example, in the fourth order approximation of V±(t) 
(see Fig. 31 the TCL as well as the NZ solution leave the 
Bloch sphere, i.e. for times larger than some critical time 
these solutions do not represent true density matrices 
anymore. 

If we look at higher orders, the NZ method is seen to be 
better than the TCL method as far as the 3-component 
of the Bloch vector is concerned. An example is shown 
in Fig.0 where we plot the tenth order approximations. 
We observe that the solution vs(t) of the TCL equation 
II58J1 is always greater than zero. This fact is obviously 
connected to the structure of this equation, which takes 
the form i?3(i) = K.3(t)v3(t) with a real function K.3(t). 
If the 3-component 1*3 of the Bloch vector vanishes at 
the time t = to, then a time-local equation of motion of 
this form can only be fulfilled if «3(<o) is also zero. In 
our case, however, the exact solution passes zero with a 
non- vanishing time derivative. It is a well-known fact& 
that the perturbation expansion of the TCL generator 
only exists, in general, for short and intermediate times 
and/or coupling strengths. This is reflected in the fact 
that the superoperator on the right-hand side of Eq. i'.Vll) 
cannot be inverted for all times, i.e. it is not always possi- 
ble to express ^3(0) in terms of v 3 (t). A similar situation 
also occurs in open systems interacting with a bosonic 
reservoir, e. g., in the damped Jaynes Cummings model 
which describes the interaction of a qubit with a bosonic 
reservoir at zero temperature. The NZ technique does 
not lead to such problems. However, since the compo- 
nents v±(t) do not vanish, the corresponding high-order 
TCL approximation is still more accurate than the NZ 
approximation, as may be seen from Fig. [3 



IV. CONCLUSION 

With the help of a simple analytically solvable model of 
a spin star system, we have discussed the performance of 
projection operator techniques for the dynamics of open 
systems and the resulting perturbation expansions of the 
equations of motion. The model consists of a central spin 
surrounded by a bath of spins interacting with the cen- 
tral spin through a Heisenberg XX coupling, and shows 
complete relaxation and partial decoherence in the limit 
of an infinite number of bath spins. Due to its high 
symmetry the model allows a direct comparison of the 
Nakajima-Zwanzig (NZ) and of the time-convolutionless 
(TCL) projection operator methods with the exact solu- 
tion in analytical terms. 

While the Born-Markov limit of the equation of mo- 
tion does not exist in the model, the dynamics of the 
central spin exhibits a pronounced non-Markovian be- 
havior. It has been demonstrated that both the NZ and 
the TCL techniques provide good approximations of the 
short-time dynamics. In practical applications the TCL 
method is usually to be preferred since it leads to time- 
local equations of motion in any desired order with a 
much easier mathematical structure, whose integration 




Figure 7: The TCL and the NZ approximation of the com- 
ponents of the Bloch vector in tenth order for a bath of 100 
spins. 



is much simpler than that of the non-local equations of 
the NZ technique. 

It should be kept in mind, however, that the expansion 
based on the TCL method converges, in general, only for 
short and intermediate interaction times. For large times 
the perturbation expansion may break down, which has 
been illustrated in our model to be connected to zeros of 
the components of the Bloch vector. It turns out that 
the NZ equation of motion yields a better approximation 
of the exact dynamics in this regime. 

In view of the heuristic approach to the Born and to 
the Redfield equation (see Sec. IIII All it is sometimes 
conjectured that a non-local equation of motion should 
be generally better than a time-local one. The results 
of Sec. IIII Bl show that this conjecture is, in general, 
not true. The fact that in a given order the time-local 
TCL equation is at least as good (and much simpler to 
deal with), if not better than the non-local NZ equation 
has also been observed in other specific system-reservoir 
models^, and has been confirmed by general mathemat- 
ical arguments^. However, care must be taken when 
applying a certain projection operator method to a spe- 
cific model: The quality of the corresponding perturba- 
tion expansion of the equation of motion may strongly 
depend on the specific properties of the model, e. g., the 
interaction Hamiltonian, the interaction time, the envi- 
ronmental state and the spectral density. 
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For example, in our particular model the TCL expan- 
sion to fourth order turns out to be more accurate than 
the fourth order NZ expansion. However, there is no rea- 
son why TCL should be generally better then NZ. To 
clarify further this point we consider the Taylor series of 
the 3-component of the Bloch vector: 



v 3 (t) =a + a,2 {at) 2 + a A (at) 4 + O ({at) 6 ) 



(61) 



The corresponding expansion obtained from TCL2 is 
given by 

a 2 

v 3 {t) = a Q + a 2 {at) 2 + -±- (erf) 4 + O ({at) 6 ) , (62) 

while NZ2 gives the expansion 

v 3 (t) = a + a 2 {at) 2 + p- {at) 4 + O ({at) 6 ) . (63) 
oao 

In our model the exact coefficients of the expansion l(6T)l 
are found to be 

16 

ao = l, a 2 = -4N, a 4 = —N 2 . (64) 

Of course, the second order coefficient a 2 is the same in 
all expansions, while in general neither TCL2 nor NZ2 
reproduce correctly the fourth order coefficient 04. How- 
ever, in our model it turns out that the TCL2 approxima- 
tion is more accurate because the fourth order coefficient 



"2 
2a 



16 9 
—N 2 
2 



(65) 



found from the solution of the TCL equation is closer 
to the correct fourth order coefficient 04 than the corre- 
sponding coefficient 



6a 



16 

if 



N 



(66) 



of the NZ equation (see Fig. EjJ. Thus we see that it 
depends crucially on the value of 04 whether TCL2 or 
NZ2 is better. 

Choosing an appropriate interaction Hamiltonian and 
initial state, one can easily construct examples where 
NZ2 is better than TCL2. For example, if vs{t) was a 
cosine function a cos(crf), then NZ2 would already give 
the exact solution. On the other hand, if v${t) was a 
Gaussian function ao exp(— a 2 t 2 ), then TCL2 would re- 
produce the exact solution because the higher cumulants 
of a Gaussian function vanish. 

The features discussed above should be taken into ac- 
count in applications of projection operator methods to 
specific open systems. In the general case in which an 
analytical solution is not known a careful analytical or 
numerical investigation of the higher orders of the respec- 
tive expansions is thus indispensable to judge the quality 
of the TCL or the NZ method, the influence of initial cor- 
relations, or to estimate the timescale over which one can 
trust the approximation obtained within a given order. 



Appendix A: BATH CORRELATION FUNCTIONS 

In this appendix we outline how to calculate the bath 
correlation functions 

Qk = ^tr B {(J + J_) fc }, (Al) 

fl*-' = ^L t r B {(J + J_) fc - Z (J_J + ) Z }. (A2) 

The trace can be computed in the eigenbasis of J3 and 
J 2 (see Sec. Ill Bll yielding a sum of polynomials in j and 
to. However, it turns out that it is easier to use the 
computational basis of the spin bath consisting of the 
states 

|si)®|s2)<8---<8> \s N ), (A3) 
where the Sj take on the values or 1 and 



4 ] \si) = (-l) Sl \Si) 



(A4) 



With the help of these states the problem is reduced to 
a combinatorial one. Since 



JV 



± 



(A5) 



we have 



{J+J-f = 4 l)(7 - 2)(J + 3)f7 - l4) • • • C7 (i2fc - l) - (i2 * ) 

il,--;i2k 

(A6) 

where the summation is taken over all possible combina- 
tions of the indices ii,i 2 , . . . ,i 2 k- Under the trace over 
the bath we can sort these indices, without interchanging 
the operators belonging to the same index, and calculate 
the partial traces over the various bath spins separately. 

Let us denote the partial trace over the Hilbert space 
of the ith bath spin by tr^. For example, we have for 
k = 2: 



tr^VM^} 
= trxjaW^ltrsI^}^!^ } 



,N-3 



0. 



(A7) 



since trj {^i^} = 0- Note that the factor 2 N 3 appears 

due to (N — 3) factors of trj {/} = 2. These factors arise 
from the partial traces of the unit matrices I in the spin 
spaces, which we did not write explicitly. As a further 
example, we have for k = 4: 

+ f (1) (3) (1) (1) (3) (1)1 

trs i <J + o_ a\' a_ a\' o_ > 

4, f (1) (1) (1) 4. f (3) (3)\. 

= tri < g + 'g + 'g_'g_' > tr3 < g_ g + ' > . 



,jV-2 



0, 



(A8) 
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because of (T^)cr±) = 0. An example of a non- vanishing 
term is given by 

, r (i) (2) (2) (i)\ 



oJV-2 



,N-2 



(A9) 



where we have used that tr* |o"^<t±H = 1. 

In view of these considerations we are now left with the 
combinatorial problem of determining all nonzero sum- 
mands for the given values of k and I. As an example, let 
us calculate explicitly the correlation function Q 2 ■ From 
its definition we have 



Q-2 



1 



N tT B {J+J-J+J-} 



1, E tr B {a^*^a^*^}.(AlO) 



It should be clear that the above method of determin- 
ing the correlation functions is easily translated into a 
numerical code from which one obtains the Qk and the 



R, 1 in any desired order. 



For N > k, the term of leading order in N of the poly- 
nomials Qk and R k ~ l is represented by the summands 
with a maximal number of k different indices, because 
these terms have the largest combinatorial weight. Af- 
ter sorting the spin operators, these terms will have the 
following form, 



■a [ i k) a { i h) . 



(A12) 



The nonzero summands in this expression have the fol- 
lowing structure: 



(*) M (*) (*) 

(J_j_ (J_ (7_j_ (T_ 
0') CO U) U) 

(*) U) U) (i) 
a\ cr_ ctx a_ 



N possibilities, 

N(N — 1) possibilities, 

N(N - 1) possibilities, 



where i ^ j in the second and the third line. Collecting 
these results we find 



There are ( k ) different ways of assigning the indices 
«i, iii ■ ■ ■ > *fc to this term. For a fixed set of indices, there 
are k\-kl different terms in the sum l|A6|l which lead to the 
sorted expression l|A12JI . corresponding to a permutation 
of the labels of all <r+ operators and of all er_ operators. 
The trace of the expression ||A12| yields 2 N ~ k . Hence, 
the term of leading order of the polynomial Qk is found 
to be 



h = ^w(N-2 N - 1 + 2-N(N-l)-2 N - 2 ) 



(All) 



-N 



N 



k\ ■ k\2 N ~ k 



N k kl 



(A13) 



A similar procedure must be carried out to calculate 



R 



k-l 



We state some results: 



= —N TV 



-iV d 



A similar proof holds for R l 1 . Thus, we have for N — > oo 
and k fixed: 



-2N + 5N A 



4N J 



-TV 4 , 
2 ' 



R 2 



-N - 



-N 



-N , 
2 ' 

4 4 ' 
19 



iiv + -iv 2 -— iV 

2 4 4 



3 , l N i 
2 ' 



Qk ~ R 



k-l 



M 
2 1 " 



N 



which has been used in Sec. Ill CI 



(A14) 
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